ELSEVIER 


Available  online  at  www.sciencedirect.com 

ScienceDirect 

Journal  of  Power  Sources  175  (2008)  712-723 

www.elsevier.com /locate /jpowsour 


JOURNAL  OF 


Microstructure  changes  in  the  catalyst  layers  of 
PEM  fuel  cells  induced  by  load  cycling 
Part  II.  Simulation  and  understanding 

Feng  Rong,  Cheng  Huang*,  Zhong-Sheng  Liu,  Datong  Song,  Qianpu  Wang 

Institute  for  Fuel  Cell  Innovation,  National  Research  Council  Canada,  4250  Wesbrook  Mall, 

Vancouver,  BC,  Canada  V6T 1W5 

Received  13  July  2007;  received  in  revised  form  1  October  2007;  accepted  1  October  2007 
Available  online  9  October  2007 


Abstract 

The  micro  structure  of  catalyst  layers  (CLs)  is  a  naturally  random  medium  and  changes  in  it  greatly  affect  the  performance  of  polymer  electrolyte 
membrane  fuel  cells.  In  this  paper,  the  mechanical  analysis  method,  developed  in  Part  I  for  understanding  the  mechanism  of  micro  structure  changes,  is 
further  extended  to  describe  CLs  as  random  three-phase  microstructures.  Microstructure  reconstruction  is  accomplished  using  statistical  information 
from  experimental  images  of  practical  CLs.  In  the  microscopically  complex  reconstructed  microstructure,  mechanical  analysis  is  performed  in 
order  to  understand  the  mechanism  of  changes  caused  by  the  cycling  of  start-up  and  shutdown  during  operation.  Numerical  simulation  shows 
that,  although  different  reconstructed  microstructures  have  different  changes,  there  have  in  common  the  competition  between  crack  initiations  in 
phases  and  delamination  between  different  phases  in  the  CLs.  This  competition  plays  an  important  role  in  microstructure  change  and  results  in 
performance  degradation,  indicated  by  the  decrease  in  connection  length  among  different  solid  components  in  the  CLs  after  certain  duty  cycles. 
Crown  Copyright  ©  2007  Published  by  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  lifetime  of  the  components  of  polymer  electrolyte 
membrane  fuel  cells  (PEMFCs)  has  a  significant  impact 
on  the  commercial  viability  of  PEMFCs  for  both  station¬ 
ary  and  transportation  energy  applications  [1].  Improvement 
in  their  durability  can  effectively  reduce  the  cost  of  imple¬ 
menting  PEMFC  systems.  Recently,  experimental  studies  have 
been  undertaken  to  understand  the  degradation  mechanism 
of  PEMFCs.  Many  researchers  [1-4]  believe  that  microstruc¬ 
ture  changes  in  catalyst  layers  (CLs)  can  be  factors  in  fuel 
cell  performance  reduction.  According  to  experiment  results, 
micro  structure  changes  in  CLs  include  the  following:  cracks, 
loss  of  carbon- supported  catalyst  clusters,  dissolution  of  recast 
electrolyte  (Nafion  ionomer),  catalyst  particle  migration,  and 
agglomerate  coarsening  [5]  These  phenomena  are  obvious 
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especially  under  dynamic  operating  conditions  (transportation 
applications).  So  the  microstructure,  especially  that  of  CLs, 
is  critical  to  durability  improvement  of  PEMFCs  [4].  There 
is  an  urgent  need  to  understand  the  underlying  mechanism  of 
microstructure  changes  in  CLs  and  how  their  evolution  reduces 
performance. 

Microstructure  changes  may  occur  in  several  ways,  such  as 
chemical  degradation  of  the  ionic  conducting  parts  or  mechan¬ 
ical  failure  in  CLs  [5].  Although  chemical  degradation  is  a  key 
factor  in  micro  structure  changes,  it  has  been  suggested  that 
mechanical  damage  is  also  highly  important  [6].  In  our  com¬ 
panion  paper  [7],  a  mechanical  analysis  model  is  developed, 
based  on  the  finite  element  method,  to  investigate  microstruc¬ 
ture  changes  in  CLs.  It  is  found  that  delaminations  and  cracks 
occur  between  phases  in  CLs  as  the  effect  of  duty  cycles.  The 
simulation  presented  in  our  companion  paper  [7]  is  based  on 
a  representation  of  microstructure  in  CLs,  which,  although  it 
offers  a  basic  understanding  of  the  changes,  suffers  from  many 
morphological  and  associated  physical  limitations.  In  addition 
to  the  simplicity  of  the  structure,  statistical  information  may 
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deviate  from  the  change  in  a  realistic  random  porous  CL.  For 
instance,  the  connection  area  between  phases  needs  to  be  fur¬ 
ther  enhanced  and  the  components  are  not  as  complex  as  those 
in  a  practical  three-phase  medium  [8-10].  To  better  understand 
the  mechanism  of  practical  random  micro  structure  changes,  a 
reconstruction  of  random  micro  structure  should  be  utilized  sys¬ 
tematically. 

In  fact,  the  reconstruction  of  random  media  from  limited 
morphological  information  (statistical  distributions,  etc.)  is  an 
intriguing  inverse  problem.  It  is  useful  to  reconstruct  CLs  using 
information  obtained  from  a  two-dimensional  (2D)  micrograph 
or  image.  A  number  of  approaches  have  been  taken  to  reconstruct 
random  media  [11,12].  Most  of  them  use  the  filtering  method 
[13]  or  stochastic  optimization  [14].  To  take  full  advantage  of 
both  of  these  methods,  we  will  apply  them  together  to  recon¬ 
struct  the  microstructure  based  on  the  statistical  information  of 
CL  images. 

First,  the  reconstruction  procedure  is  proposed:  statistical 
feature  extraction  based  on  an  image  process  of  CLs,  initial 
reconstruction  based  on  the  filtering  method,  and  stochastic  opti¬ 
mization  to  refine  the  results  of  the  filtering  method.  Second, 
the  method  developed  in  our  companion  work  [7]  is  extended  to 
investigate  the  reconstructed  micro  structure  changes.  Finally, 
a  basic  understanding  of  the  mechanism  of  microstructure 
changes  in  CLs  is  achieved  through  mechanical  analysis.  More¬ 
over,  a  common  indication  of  performance  degradation  for 
different  reconstruction  results  is  observed  in  the  simulation 
results. 

2.  Microstructure  reconstruction  in  CLs 

The  determination  of  the  macroscopic  behaviour  of  CLs,  a 
random  porous  medium,  can  be  performed  in  two  steps:  gen¬ 
erating  media  with  specific  geometrical  or  statistical  properties 
and  solving  the  relevant  set  of  equations.  A  reconstruction  pro¬ 
cess  based  on  2D  transmission  electron  microscopy  (TEM) 
experimental  images  is  proposed.  In  this  paper,  the  whole  recon¬ 
struction  of  the  microstructure  can  be  divided  into  three  stages: 
evaluation  of  phase  statistical  information,  initial  reconstruc¬ 
tion  based  on  the  filtering  method,  and  stochastic  optimization 
reconstruction  of  the  microstructure. 

2.7.  Evaluation  of  phase  statistical  information 

A  common  starting  point  of  the  reconstruction  procedure  is 
the  statistical  distribution  and  correlation  relationship  of  phases. 
For  a  CL,  the  phases  include  the  electrolyte  (Nation),  the  C/Pt 
agglomerate  and  pores.  We  want  to  generate  a  random  porous  CL 
with  volume  fractions  and  correlation  functions  of  components. 
The  CL  is  assumed  to  be  homogenous  and  isotropic,  although  the 
latter  is  not  essential.  First,  experimental  images  are  processed  to 
get  the  statistical  information.  In  this  section,  one  method,  based 
on  the  physical  properties  of  the  different  phases  in  CLs,  is  pro¬ 
posed  to  segment  the  phases  in  2D  experiment  results.  In  black 
and  white  figures,  this  information  is  conveyed  by  the  gray  value 
or  the  intensity  of  different  pixels.  So  we  assume  that  every  pixel 
represents  an  individual  phase  and  the  intensity  of  the  pixel  can 


be  used  to  indicate  its  properties.  Then  three  intervals  of  inten¬ 
sity  from  0  to  255  should  be  determined  to  distinguish  different 
phases. 

The  intensity  of  Nafion  is  assumed  to  be  greater  than  that 
of  the  C/Pt  agglomerate  and  lower  than  that  of  the  pores.  To 
determine  these  two  boundary  values,  the  physical  properties 
and  mass  ratios  of  the  different  phases  must  be  used.  When  the 
TEM  or  scanning  electron  microscopy  (SEM)  images  of  a  CL  are 
taken,  the  mass  ratio  £,  as  well  as  the  density  p,  of  the  different 
phases  is  known.  So  the  volume  ratio  (8)  of  the  different  phases 
can  be  determined  as  follows: 


8\  =  1  —  82  —  <$3,  82 


<$3 


PCL^C  PClAPt 
PC  PPt 


PCL^Nafion 

PNafion 


(1) 


where  82,  ^Nafion,  and  PNafion  are  volume  ratio,  mass  ratio,  and 
density  of  Nafion,  respectively;  sq  and  pc  are  the  mass  ratio  and 
density  of  carbon,  respectively  ;£pt  and  ppt  are  the  mass  ratio  and 
density  of  platinum,  respectively;  8\  and  <$3  are  the  volume  ratio 
of  the  pores  and  C/Pt  agglomerate,  respectively;  and  pcl  is  the 
density  of  the  whole  CL. 

When  the  CL  is  assumed  to  be  homogenous,  statistical  aver¬ 
ages  can  be  replaced  by  volume  averages.  When  it  is  assumed 
to  be  isotropic,  these  volume  averages  can  be  replaced  by  sur¬ 
face  averages.  Hence,  the  use  of  thin  sections  is  justified.  So 
the  volume  ratios,  which  can  be  regarded  as  area  ratios  in  the 
2D  experimental  figures,  can  be  used  to  evaluate  the  bound¬ 
aries  of  intensity.  Fig.  1(a),  a  typical  image  of  a  CL  by  Xie 
et  al.  [15],  is  used  to  illustrate  the  process  of  determining 
boundaries  of  intensity.  In  their  fabrication,  catalyst  inks  were 
prepared  by  mixing  catalyst  powder  (20wt%  Pt  on  a  Vulcan 
XC-72,  E-TEK),  Nafion  solution,  and  isopropanol.  The  CL  in 
Fig.  1(a)  contains  30wt%  Nafion  115.  In  fact,  when  Xie  et 
al.  [15]  were  verifying  the  performance  with  the  macrohomo- 
geneous  model,  the  volume  fraction  of  different  compositions 
in  the  CL  were  reproduced  based  on  the  constitutive  relation¬ 
ship  between  weight  fractions  and  volume  fractions  (Eq.  (1)); 
that  is,  C/Pt:Nafion:pore  =  46.7:25.7:27.6.  The  distribution  of 
the  intensity  of  different  pixels  is  shown  in  Fig.  1(b).  7a  and  7b 
are  the  segment  phases  in  the  image.  Based  on  the  assumption 
of  isotropic  random  media,  area  fractions  of  the  CL  components 
should  be  the  same  as  the  volume  fractions  in  the  experiments. 
So  7a  and  7b  are  searched  for  from  0  to  255  to  best  approach 
these  area  fractions.  7a  is  determined  to  be  94  and  7b  is  119.  In 
fact,  the  area  fractions  of  the  components,  based  on  intensity,  are 
43.4:26.7:30.3  (C/Pt:Nafion:pore).  When  these  two  values  have 
been  determined,  Fig.  1(a)  can  be  filtered  to  a  black-and-white 
figure  to  indicate  different  phases,  as  shown  in  Fig.  2  where  white 
pixels  are  the  reference  phases,  specifically  (a)  pore,  (b)  Nafion, 
and  (c)  C/Pt  agglomerate.  At  this  point,  the  image  segmentation 
has  been  accomplished  based  on  the  mass  ratios,  densities  of 
components,  and  intensity  of  pixels. 

As  mentioned  previously,  the  presence  of  different  phases  in 
the  sample  is  measured.  At  each  point  (pixel)  within  the  experi¬ 
mental  sample,  an  overall  phase  function  Z(x)  at  each  point  x  is 
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100  nm 

Fig.  1.  (a)  The  TEM  image  of  a  cathode  CL  from  Ref.  [15]  and  (b)  intensity  distribution  of  the  TEM  image. 


250 


Fig.  2.  The  segmentation  of  the  TEM  image  (Fig.  1(a))  based  on  intensity,  where  the  white  pixels  are  (a)  pores,  (b)  Nation,  and  (c)  C/Pt  agglomerates. 
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defined  as 


Z(x)  =  < 


co  l 
002 
(02 


if  x  belongs  to  pore 

if  x  belongs  to  Nation 

if  x  belongs  to  C/Pt  agglomerate 


(2) 


The  set  of  values  { (Ok ,  k=  1,  2,  3}  is  arbitrary;  however,  in 
this  paper,  this  set  is  equal  to  the  sequence  of  the  first  three  inte¬ 
gers  {0,  1,  2}.  As  x  varies,  Z(x)  describes  a  discrete  stochastic 
process.  Assuming  a  homogenous  and  isotropic  medium  implies 
that  this  process  is  stationary;  that  is,  it  has  a  constant  mean  and 
its  autocovariance  function  is  invariant  under  arbitrary  transla¬ 
tions.  Theoretically,  the  overall  average  of  Z(x)  and  the  two-point 
overall  correlation  function  can  be  introduced  as 


X  =  Z(x) 


Rz(  u)  = 


(Z(x)  ~  x)(Z(x  +  u)  -  x) 
(Z(x)2  -  x2) 


(3) 


where  Rz(u)  is  the  two-point  correlation  function,  which  rep¬ 
resents  the  probability  that  two  points  at  a  distance  x  are  both 
the  same  phase;  the  overbar  denotes  the  statistical  average;  and 
u  is  the  lag  vector.  For  a  statistically  homogenous  CL,  x  is  a 
constant  and  Rz(u)  is  a  function  only  of  the  lag  vector,  u;  that 
is,  it  is  independent  of  the  location  vector  x.  Furthermore,  if  the 
CL  is  isotropic,  the  correlation  function  Rzivf)  is  a  function  only 
of  the  norm  of  u. 

This  is  a  theoretical  declaration  of  the  overall  mean  and  two- 
point  correlation  function.  It  is  still  necessary  to  mention  the 
procedure  by  which  they  are  computed  from  2D  experimental 
images  of  an  actual  CL.  To  minimize  finite  size  effects,  periodic 
boundaries  and  more  experimental  figures  are  utilized  during 
evaluation.  In  Fig.  2,  the  M  x  N  ( 781  x781)  pixels  2D  image  is 
defined  as  a  discrete  valued  function  Z(v,  y),  where  Z(x ,  y)  is  the 
same  as  Z(x)  in  Eq.  (2).  The  three-phase  image  S  is  divided  into 
two  halves,  S i  and  S2,  which  satisfy  S1US2  =S  and  S\HS2  =  0. 
To  calculate  /?z(u),  S 1  is  first  translated  by  a  distance  u  along 
the  x-axis;  it  yields  S\(+u).  The  spatial  average  indicated  in  Eq. 
(3)  is  replaced  by  an  intersection  of  images: 


Z(x,  y)Z(x  +  u,y)  =  S\(+u)  (T  S 


The  other  operations  indicated  in  Eq.  (3)  are  then  performed 
algebraically.  Fig.  3  shows  the  result  from  the  three-phase 
images  in  Fig.  2.  General  practical  recommendations  for  the 
measurements  of  these  correlation  functions  can  be  found  in 
Ref.  [13]. 

In  addition,  it  is  possible  to  calculate  numerically  the  inter- 
and  auto-correlations  between  the  different  phases  from  the 
three-phase  images.  Theoretically,  the  three-phase  functions 
Zk(x)  can  be  defined  as 

f  1  if  Z(x)  =  (ok 
Zk(x)  =  < 

I  0  otherwise 

The  expected  value  (denoted  by  an  overbar)  of  Z&(x)  is  the 
probability  of  the  presence  of  this  phase: 

Xk  =  Z*(x) 


Fig.  3.  The  global  correlation  function  of  the  TEM  image  (Fig.  1(a)). 


Here  Xk  should  be  equal  to  8k  in  Eq.  (1)  if  7a  and  4  are 
determined  from  8k.  Notice  that  x  =  J2k=l ^kXk-  The  correla¬ 
tion  function  Rzk,zm(}x)  is  the  two-point  correlation  function, 
which  is  the  probability  that  one  point  at  a  distance  x  belongs 
to  phase  k  and  one  point  at  a  distance  x  +  u  is  phase  m.  The 
definition  of  Rzk,zm(u 0  is 

(Zjt(x)  -  X£)(Zm(x  +  u)  —  Xrn) 

Rzk,zjn)  = - -  — - 

y(  Xk  +  XjfcXXm  +  xi) 

(k,m  =  1,2,  3)  (4) 

So  from  three-phase  images  (Fig.  2),  the  cross-  and  auto¬ 
correlation  relationship  can  be  calculated,  as  shown  in  Fig.  4 
where  (a)  is  the  auto-correlation  functions  and  (b)  is  the  cross¬ 
correlation  functions  between  phases. 

In  this  section,  the  segmentation  of  a  gray  image  of  an  actual 
CL  is  performed  based  on  mass  ratios  and  physical  constants  of 
the  components.  In  addition,  the  global  correlation  function,  as 
well  as  functions  between  phases,  can  be  calculated. 

2.2.  Initial  reconstruction  based  on  the  filtering  method 

The  reconstruction  of  micro  structure  in  CLs  can  be  achieved 
based  on  feature  extraction  of  an  actual  three-phase  CL.  We  seek 
to  generate  a  random  porous  microstructure  containing  three 
phases,  each  having  a  specified  volume  fraction  8k.  The  over¬ 
all  correlation  function  Rz(u)  is  given  (Eq.  (3)  and  Fig.  3).  An 
extensively  examined  reconstruction  method,  called  the  filtering 
method,  is  based  on  successively  passing  a  normalized  uncor¬ 
related  random  Gaussian  field  through  first  a  linear  and  then 
a  nonlinear  filter  to  yield  the  discrete  values,  which  represent 
the  phases  of  the  structure.  In  this  section,  the  filtering  method 
is  used  to  reconstruct  an  initial  value  for  subsequent  stochastic 
optimization.  The  reason  for  this  will  be  discussed  at  the  end  of 
this  section. 

Reconstructing  a  three-phase  medium  is  done  by  thresh¬ 
olding  a  continuous  Gaussian  field.  The  initial  reconstructed 
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Fig.  4.  The  cross-  and  auto-correlation  relationships  of  TEM  image  (Fig.  1(a)). 


microstructure  consisting  of  a  Gaussian  distributed  random  field 
X  is  generated  using  Marsaglia’s  ziggurat  algorithm  [16],  as 
shown  in  Fig.  5(a).  This  independent  Gaussian  field  X  is  con¬ 
voluted  by  one  constructed  linear  filter  and  forms  another  field 
Y  that  is  still  Gaussian  distributed  but  correlated.  The  nonlinear 
filter  then  performs  a  threshold  cut  to  the  field  Y  to  generate  the 
final  reconstructed  structure  Z.  The  following  will  gives  the  con¬ 
struction  and  influence  of  these  filters  and  relates  their  properties 
to  the  statistical  properties  of  the  resulting  fields.  The  construc¬ 
tion  of  filters  is  called  an  inverse  problem  by  Losic  et  al.  [17]. 
For  practical  purposes,  the  porous  microstructure  is  constructed 
in  a  discrete  manner.  It  is  considered  to  be  composed  of  Nx  x  Ny 
(in  Fig.  5(a),  Nx  =  Ny  =  200)  small  squares,  each  of  the  same  size 
a.  These  elementary  squares  are  filled  with  either  voids,  Nafion, 
or  C/Pt  agglomerate.  Hence,  the  spatial  variables  x  and  u  in  Eqs. 
(3)  and  (4)  will  take  only  discrete  values;  the  corresponding  trios 
of  integers  are  denoted  by 

x  =  (/,  j)  *  a  and  u  =  (r,  s)  *  a ,  i,  re[  1,  Nx]; 
j,se[l,Ny]  (5) 

Moreover,  a  set  of  strictly  increasing  constants  {£>;,  z  =  0,  1, 
. . .,  n}  with  bo  =  0  and  bn  =  1  is  constructed.  It  partitions  the 
interval  [0,  1]  into  n  segments  Ii  =  [bi^%,  bi\ ,  7=1,  2,  . . .,  n. 
Accordingly,  the  set  91  of  real  numbers  is  partitioned  into  the 


family  of  interval 

Ji  =  [0-\bi-l),0-\bi)]  (6) 

where  <P(y)  is  the  Gaussian  probability  distribution  function 
(PDF);  that  is, 


A  value  of  cok  among  the  possible  values  of  Z  in  Eq.  (2)  is 
associated  with  each  of  these  segments  //;  that  is,  each  part  is 
considered  to  belong  to  one  of  three  phases.  Note  that  n  can  be 
larger  than  three,  but  is  not  necessarily  equal;  several  zones  can 
correspond  to  the  same  phase. 

The  most  difficult  point  is  the  determination  of  the  linear  fil¬ 
ter.  One  can  start  from  the  fact  that  the  random  vector  (F(x), 
Y(x  +  u)),  after  a  filter  is  applied,  is  a  bivariate  Gaussian  whose 
probability  density  is  known;  this  density  can  be  expanded  in 
terms  of  Hermite  polynomials.  After  some  systematic  manipu¬ 
lations,  which  use  classical  identities  [18],  the  target  correlation 
function  Rz(vl)  can  be  expressed  as  a  series  in  terms  of  the 
correlation  function  of  Y,  Ry(u): 

00 

Rz(  U)  =  J2c2mRy( u)  (8) 

m= 0 


Fig.  5.  (a)  Uncorrelated  Gaussian  distributed  random  field.  Here  Nx  =  Ny  =  200;  (b)  coefficients  of  linear  filter;  (c)  correlated  Gaussian  distributed  random  field  after 
linear  filter  applied. 
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The  coefficients  Cm  are  given  by 


greatly  reduced  since  c  is  now  a  function  of  only  the  distance 


1  “  r<P~l(bi) 

- —  /  ( 

•7tm\crz  1) 


Cn 


V2 


i~(y2^Hm(y)dy 


where  <P(y)  is  the  Gaussian  PDF  (Eq.  (7)),  bi  is  the  strictly 
increasing  constants  set  defined  in  Eq.  (6),  az  is  the  square  root 
of  the  variance  of  the  random  variable  Z  and  Hm(y )  are  Hermite 
polynomials.  They  may  be  expressed  as 

°\  =  (y,n<7\  -  ( 


».W  =  (-ir«P  (y)  («P  (-y)) 

When  Xk  is  given  from  Section  2.1,  the  correlation  function 
Ry( u)  is  easily  derived  from  7?z(u),  based  on  Eq.  (3).  This  simply 
corresponds  to  the  numerical  inversion  of  Eq.  (8)  by  any  stan¬ 
dard  method.  In  our  reconstruction,  an  algorithm  developed  by 
Dekker  [19]  uses  a  combination  of  bisection,  secant,  and  inverse 
quadratic  interpolation  methods. 

Once  Ry( u)  is  known,  the  linear  filter  may  be  constructed.  A 
linear  operator  can  be  defined  by  an  array  of  coefficients  c(v), 
where  v  belongs  to  a  finite  square  [0,  Lc ]2  in  Z2.  Outside  this 
square,  it  is  equal  to  0.  A  new  random  field  7(x)  can  be  expressed 
as  a  linear  combination  of  the  uncorrelated  Gaussian  random 
variables  X(x), 

Y(x)=  J2  c(v)X(x  +  v)  (9) 

v  e  [0,LC]2 

Their  correlation  function  Ry( u)  is  easily  seen  to  be 
fly(u)  =  Y,  c(v)c(v  + u)  (10) 

v  e  [0,LC]2 

if  the  variance  of  Y(x)  is  equal  to  1 . 

For  isotropic  media  and  the  discretization  of  microstructure 
(Eq.  (5)),  the  number  of  coefficients  c  in  the  linear  filter  can  be 


c(v)  =  c(r,  s)  =  c(\/ r2  +  s2)  =  c(d)  (11) 

where  d  is  the  distance.  Substituting  Eq.  (11)  into  Eq.  (10),  one 
can  determine  the  function  c  from  nonlinear  equations: 

Y  c(\/r2  +  s2)c(\J  (r  +  u)2  +  s2)  =  Ry{u)  (12) 

ve[0,Lc]2 

In  our  simulation,  the  nonlinear  least  squares  curve-fitting 
algorithm  is  utilized  to  obtain  the  coefficients  c  in  the  linear 
filter  from  Eq.  (12).  The  coefficients  c  are  shown  in  Fig.  5(b). 

The  random  field  y(x)  (shown  in  Fig.  5  (c))  is  correlated  after 
the  linear  filter  (Eq.  (9))  is  applied.  But  this  is  still  not  satisfactory 
since  it  takes  its  value  in  91,  while  the  CL  must  be  represented 
by  a  three- valued  field  Z(x).  In  order  to  extract  such  a  field  from 
F(x),  a  nonlinear  filter  G  is  applied;  that  is,  the  random  variable 
Zis  a  deterministic  function  of  Y,  Z=G(Y).  When  the  value  y  of 
the  random  variable  Y(x)  falls  within  the  interval  7/  (Eq.  (6)),  Z 
is  set  to  the  corresponding  value  cot 

Z  =  G(Y)  =  (Of  if  ye  Jt 

Obviously,  the  probability  of  each  coi  is  equal  to  Xk  if  the 
following  condition  is  fulfilled 

yf  (bi  -  bi- 1)  =  Xk 

o)i=cok 

At  this  point,  the  nonlinear  filter  is  constructed  to  correlate 
the  random  field  in  91  to  the  three-phase  field  Z,  which  is  the 
reconstruction  of  the  micro  structure.  Fig.  6(a)  shows  the  result 
after  the  nonlinear  filter  is  applied  to  Fig.  5(c).  The  overall  corre¬ 
lation  function  of  Fig.  6(a)  can  be  evaluated  based  on  the  method 
proposed  in  Section  2.1,  as  shown  in  Fig.  6(b).  It  can  be  seen 
that  the  difference  between  the  reconstructed  and  target  corre¬ 
lation  functions  is  relatively  small,  especially  in  the  short-range 
correlation. 

In  the  previous  filtering  method,  the  only  correlation  imposed 
on  the  process  is  the  global  correlation  function  (Eq.  (3)).  Hence, 
the  cross-correlation  relationship  between  the  different  phases 
in  the  reconstructed  media  should  be  examined  to  verify  the 


Fig.  6.  (a)  Three-phase  reconstructed  micro  structure  (black  pixels  represents  C/Pt  agglomerate,  gray  pixels  are  Nation  and  white  are  pore);  (b)  overall  correlation 
function  of  reconstructed  microstructure. 
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Fig.  7.  The  auto-correlation  relationship  of  pore  in  both  reconstructed 
microstructure  and  experimental  image. 


the  calculated  and  target  functions  via  a  stochastic  optimization 
technique,  such  as  the  simulated  annealing  method  [14]. 

Although  the  filtering  method  has  several  disadvantages,  it 
gives  us  a  good  starting  point  for  reconstruction.  So  in  this  paper, 
stochastic  optimization  starts  from  the  reconstruction  of  the  fil¬ 
tering  method.  The  advantage  of  this  initial  value  of  optimization 
will  be  presented  at  the  end  of  this  section.  Now,  we  will  focus  on 
the  key  points  of  stochastic  optimization  for  the  reconstruction 
of  microstructures  in  CLs. 

Consider  the  reference  (experimental)  three-phase  cross¬ 
correlation  functions  Rzk,zm  (Km  =1,2,  3)  of  a  CL  presented  in 
Section  2. 1  and  calculated  from  the  2D  experimental  images.  Let 
Rzk,zm  be  the  corresponding  correlation  function  of  the  recon¬ 
structed  microstructure  (with  periodic  boundary  conditions)  at 
some  time  step.  It  is  this  micro  structure  that  we  will  attempt 
to  evolve  towards  Rzk,zm  from  the  reconstruction  result  of  the 
filtering  method. 

At  any  particular  time  step,  a  fictitious  energy  E  can  be  defined 
as 


results.  The  auto-correlation  of  pore  phases  is  calculated  accord¬ 
ing  to  Eq.  (4)  and  Fig.  6(a)  and  shown  in  Fig.  7.  Although 
the  trend  of  auto-correlation  functions  between  different  phases 
follows  that  of  experimental  ones,  there  are  clearly  noticeable 
discrepancies  between  them.  This  can  be  considered  to  be  one  of 
drawbacks  of  using  the  filtering  method  to  reconstruct  multiple- 
phase  microstructure. 

In  fact,  it  should  be  emphasized  that  there  are  several  arti¬ 
ficial  parameters  introduced  in  the  filtering  method:  the  size  of 
the  generated  microstructure  Nx  and  Ny,  the  cutoff  maximum 
m  in  the  infinity  series  of  the  linear  filter  (Eq.  (8)),  the  length 
of  the  linear  filter’s  square  Lc  (Eq.  (9)),  and  number  of  seg¬ 
ment  intervals  in  [0,  1]  n  for  constructing  the  nonlinear  filter. 
These  parameters  dramatically  influence  the  reconstruction  of 
the  microstructure  [13]  and  are  also  a  key  problem  of  the  filtering 
method  for  reconstruction. 

Moreover,  it  is  difficult  to  extend  the  Gaussian  filter¬ 
ing  method  to  incorporate  non-Gaussian  statistics.  Hence  the 
method  is  model-dependent;  that  is,  it  depends  on  the  underly¬ 
ing  Gaussian  statistics.  So  in  our  proposed  reconstruction,  the 
results  of  the  filtering  method  are  used  as  the  initial  value  for 
stochastic  optimization  for  reconstruction. 

2.3.  Stochastic  optimization  for  reconstruction 

Clearly,  the  filtering  method  with  the  conventional  overall 
correlation  function  alone  may  not  be  adequate  to  characterize 
the  microstructure  of  multiphase  media  for  accurate  reconstruc¬ 
tion.  It  is  desirable  that  a  reconstruction  procedure  has  the  ability 
to  incorporate  as  much  crucial  microstructure  information  as 
possible  to  capture  the  salient  features  of  the  reference  structure. 

Recently,  Torquato  et  al.  [12,20-22]  developed  a  reconstruc¬ 
tion  procedure  based  on  a  stochastic  optimization  technique. 
Starting  from  an  initial  realization  of  the  random  medium,  the 
method  proceeds  to  find  a  realization  in  which  the  calculated 
correlation  functions  best  match  the  target  functions.  This  is 
achieved  by  minimizing  the  sum  of  squared  differences  between 


E  =  EE^zJu)  -  RZk,zm(u)]2 
u  k,m 

where  the  sum  of  u  is  over  all  discrete  values  of  the  gener¬ 
ated  microstructure  (Eq.  (5))  and  the  indexes  k  and  m  denote 
the  type  of  the  different  phases.  To  evolve  the  microstructure 
towards  Rzk,zm  (he.,  minimizing  E ),  we  interchange  the  states 
of  two  arbitrarily  selected  pixels  of  different  phases,  automati¬ 
cally  preserving  the  volume  fraction  of  both  phases.  After  the 
interchange  is  performed,  the  new  energy  E'  and  the  energy  dif¬ 
ference  A E  -E!  —  E  can  be  calculated.  This  phase  interchange 
is  then  accepted  with  some  probability  p(AE )  that  depends  on 
A E.  Here  the  probability  is  given  by  the  Boltzmann  distribution; 
that  is,  the  Metropolis  acceptance  rule  [20] 


p(AE)  = 


1  AE  <  0 

exp(-A£/T)  AE  >  0 


where  T  is  a  fictitious  temperature.  This  is  the  concept  of  a  sim¬ 
ulated  annealing  schedule.  The  cooling  or  annealing  schedule, 
which  governs  the  value  and  rate  of  change  of  T,  is  chosen  to 
be  sufficiently  slow  to  allow  the  system  to  evolve  to  the  desired 
state  as  quickly  as  possible  without  trapping  in  any  local  energy 
minima  (metastable  states).  At  each  annealing  step  k ,  the  sys¬ 
tem  is  allowed  to  evolve  long  enough  to  thermalize  at  T(k).  The 
temperature  is  then  lowered  according  to  a  prescribed  anneal¬ 
ing  schedule  T(k)  until  the  energy  of  the  system  approaches  its 
ground  state  value  within  an  acceptable  tolerance. 

In  our  reconstruction  procedure,  the  system  evolves  through 
the  logarithmic  annealing  schedule,  which  decreases  the  ficti¬ 
tious  temperature  to  the  ground  state  according  to  T(k)  ~  1/ln  (k). 
A  logarithmic  decrease  may  cause  very  slow  convergence.  Thus, 
for  practical  purposes,  we  adopt  the  more  popular  and  faster 
annealing  schedule  T(k)/T( 0)  =  kk,  where  the  constant  A  is  the 
annealing  rate  which  must  be  less  than  butclose  to  1  and  assumed 
before  annealing.  More  detail  can  be  found  in  Ref.  [14]. 

As  a  matter  of  fact,  the  optimization  problem  normally 
faces  the  same  challenge,  the  sensitivity  of  the  initial  value. 
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Fig.  8.  (a)  Microstructure  reconstruction  after  stochastic  optimization  (initial  pattern  is  Fig.  12(a),  black  pixels  represent  C/Pt  agglomerate,  gray  pixels  are  Nation, 
and  white  are  pore);  (b)  “temperature”  dependence  of  configuration  energy. 


The  character  of  the  energy  landscape  determines  how  closely 
one  can  numerically  approach  the  true  global  minimum.  In 
general,  bad  initial  values  of  optimization  will  result  in  both 
local  metastable  states  and  meaningless  computation  time.  To 
avoid  this  effect,  two  kinds  of  reconstruction  are  proposed. 
One  starts  from  a  random  three-phase  pattern,  while  the  other 
is  based  on  the  result  of  the  filtering  method  (Fig.  6(a)).  In 
the  latter,  the  result  of  the  filtering  method  is  chosen  as  the 
initial  value  of  stochastic  optimization  and  it  has  been  ver¬ 
ified  that  the  overall  correlation  function  (Fig.  6(b))  recurs 
the  target  one  better.  Fig.  8(a)  shows  the  result  of  stochastic 
optimization. 

Fig.  8(b)  shows  the  “temperature”  dependence  of  configura¬ 
tion  energy  E.  The  optimization  computational  time  cost  of  the 
initial  random  pattern  is  18,061.31  s  while,  if  we  start  from  the 


result  of  the  filtering  method,  the  computational  time  is  about 
one  fourth  of  that  (i.e.,  4615.22  s).  These  reconstructions  show 
the  advantages  of  using  the  results  from  the  filtering  method. 
Doing  so  makes  the  energy  less  than  the  random  initial  value 
and  saves  computational  time. 

Moreover,  the  lower  order  correlation  functions  do  not  con¬ 
tain  complete  morphological  information  and  thus  they  cannot 
uniquely  characterize  the  micro  structure,  even  if  the  global 
minimum  is  achieved.  The  generated  microstructures  do  not  nec¬ 
essarily  contain  the  same  statistical  information  beyond  what  is 
imposed,  even  though  a  naked-eye  comparison  is  many  times 
more  insensitive  to  such  differences.  So  we  reconstruct  five 
different  microstructures  (Fig.  9,  Nx=Ny  =  20)  to  understand 
the  mechanical  mechanism  of  micro  structure  changes  due  to 
environmental  change.  According  to  Fig.  9(f),  the  correlation 


Fig.  9.  (a)-(e)  Five  different  microstructure  reconstructions  used  to  simulate  the  microstructural  changing  (Nx  =  Ny  =  20,  black  pixels  represent  C/Pt  agglomerate, 
gray  pixels  are  Nafion,  and  white  are  pore);  (f)  the  auto-correlation  relationship  of  pore  in  these  reconstructed  microstructures  and  experimental  image. 
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Fig.  10.  Different  phases  in  CLs,  where  the  dashed  line  is  the  boundary  of  Nahon 
and  the  C/Pt  agglomerate,  and  solid  lines  are  the  boundaries  of  pore  and  the  other 
two  components. 

functions  of  pores  in  these  microstructure  reconstructions  are 
good  estimators  of  that  of  the  experimental  image. 

In  summary,  the  microstructure  in  CLs  is  recurred  through 
a  reconstruction  procedure.  Firstly,  the  experimental  image  is 
analyzed  and  the  correlation  functions  of  phases  are  extracted 
as  the  target  of  reconstruction.  Secondly,  the  filtering  method 
is  utilized  to  give  us  the  initial  estimation  of  reconstruction. 
Finally,  the  stochastic  optimization  is  applied  to  get  the  final 
reconstruction. 

3.  Mechanical  analysis  of  microstructure  changes 


on  the  reconstructed  CL.  The  mechanical  analysis  model,  devel¬ 
oped  for  fictitious  microstructure  changes  in  our  companion 
paper  [7],  is  extended  here  to  investigate  the  changes  in  recon¬ 
structed  microstructure  due  to  environmental  change,  such  as 
thermal  and  humidity  cycles  during  the  start-up  and  shutdown 
of  the  fuel  cell. 

As  indicated  in  Fig.  10,  there  are  three  different  phases  in  CLs: 
electrolyte  (Nafion),  pore,  and  carbon/catalyst  (Pt)  agglomerate. 
In  our  model,  two  components  (Nafion  and  the  C/Pt  agglomer¬ 
ate)  are  meshed  in  the  domain  for  the  sake  of  the  finite  element 
method.  The  pore  phase  is  considered  as  void  in  the  solid. 
To  simulate  the  mechanical  response  of  Nafion  and  the  C/Pt 
agglomerate,  the  rate-independence  constitutive  relationship  is 
applied  to  depict  the  relationship  between  the  force  and  defor¬ 
mation  in  them.  It  is  assumed  that  the  carbon  and  platinum  have  a 
linear-elastic  behaviour  and  do  not  swell  in  response  to  moisture. 
Based  on  many  mechanical  experimental  results,  the  piecewise 
linear  isotropic  elastic-plastic-failure  model  (Fig.  11(a))  is  pro¬ 
posed  to  describe  the  behaviour  of  Nafion.  The  parameters  in  the 
model  are  Young’s  modulus  Eq,  plastic  modulus  E\,  static  yield 
stress  crpo ,  and  failure  stress  as.  These  parameters  depend  on 
changes  in  moisture  and  temperature.  Two-dimensional  inter¬ 
polation  is  used  to  obtain  all  data  between  discrete  experimental 
data  points.  Young’s  modulus  Eo,  for  example,  is  shown  in 
Fig.  11(b)  and  the  following 


Eo_ 


,  rhV  /rh, 

1.9  -  -20.2  -  +57500 

RHr  /  V  RFL 


x  |  0.023(0  -0.213^1+387 


As  indicated  in  Section  2,  the  reconstruction  process  is  based 
on  the  experimental  work  of  Xie  et  al.  [15].  Based  on  the  experi¬ 
mental  TEM  images  of  CLs  after  fabrication,  the  reconstruction 
is  processed  with  the  stochastic  optimization  method. 

Once  the  microstructure  is  reconstructed  and  the  constituent 
phases  (i.e.,  carbon/catalyst  agglomerate,  pore,  and  electrolyte 
phases)  are  identified,  the  mechanical  analysis  can  be  performed 


where  RHr,  Tr,  and  Eqy  are  the  reference  values  of  humidity, 
temperature,  and  Young’s  modulus,  respectively.  The  reference 
Young’s  modulus  is  input  as  the  value  (226.2  MPa)  at  30%  RH 
and  25  °C.  The  other  parameters  can  be  found  in  our  companion 
paper  [7]. 

In  addition,  to  simulate  the  debonding  phenomena  (or  delam¬ 
ination)  between  Nafion  and  the  C/Pt  agglomerate,  the  cohesive 


Fig.  11.  (a)  The  piecewise  linear  constitutive  relationship  of  Nafion  adopted  in  the  model;  (b)  Young’s  modulus  as  functions  of  temperature  and  relative  humidity. 
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zone  model  (CZM)  is  applied  to  mimic  the  boundary  of  these 
two  phases,  indicated  as  a  dashed  line  in  Fig.  10.  The  CZM  is,  in 
fact,  a  model  that  adopts  softening  relationships  between  forces 
and  separations,  which  in  turn  introduce  a  critical  fracture  energy 
that  is  also  the  energy  required  to  break  apart  the  interfaces.  The 
parameters  in  the  CZM  can  be  determined  through  the  empirical 
relationship  between  the  boundary  and  bulk  material.  It  should 
be  mentioned  that  the  constitutive  relation  for  each  cohesive 
surface  is  taken  to  be  elastic  so  that  any  dissipation  associated 
with  separation  is  ignored.  So  the  failure  of  CZM  elements  is 
manually  considered  when  the  strain  in  the  CZM  exceeds  one 
initial  strain  strength  for  the  statistics  of  interface  failure  between 
Nation  and  the  C/Pt  agglomerate;  that  is,  for  one  CZM  element, 
if 


max 


>  co 


(13) 


the  element  is  considered  to  be  a  failed  one.  Here  8n  is  the  normal 
separation  across  the  interface,  8t  is  the  shear  separation  along 
the  interface,  8n  is  the  normal  separation  where  the  maximum 
normal  traction  is  attained  with  8t  =  0,  8t  is  the  shear  separation 
where  the  maximum  shear  traction  is  attained  at  8t  =  8t/V^‘, 
and  co  is  the  coefficient  corresponding  to  the  maximum  strain 
of  interfaces. 

When  fuel  cells  start-up  and  shutdown,  the  humidity  and  tem¬ 
perature  in  the  CLs  will  change,  and  the  Nation  will  swell  or 
shrink.  Therefore,  the  Nation  will  come  into  contact  with  the 
C/Pt  agglomerate  or  with  another  Nation  domain.  So,  the  contact 
model  is  used  on  the  boundary  between  the  pore  and  the  other 
two  components  (solid  lines  in  Fig.  10).  All  of  these  contact  pairs 
include  the  effect  of  friction,  in  which  the  coefficient  of  friction 
is  determined  through  experiment  [23].  It  should  be  noted  that 
these  contacts  have  no  stickiness;  that  is,  when  the  normal  force 
across  the  interface  is  traction,  the  contact  phenomenon  will 
disappear. 

In  summary,  the  finite  element  method,  with  the  support  of 
the  CZM  and  the  frictional  contact  model,  is  proposed  to  obtain 
the  mechanical  response  of  CLs.  The  parameters  in  the  CZM  and 
the  contact  model  can  be  determined  from  the  empirical  relation, 
as  done  in  Part  I.  With  the  finite  element  method  simulation, 
the  underlying  mechanism  of  the  changes  of  the  reconstructed 
microstructure  can  be  understood. 


4.  Mechanism  of  reconstructed  microstructure  changes 

According  to  observations,  cracks  and  delamination  can  be 
found  in  the  aging  CLs  of  PEM  fuel  cells.  The  finite  element 
method,  along  with  the  CZM  and  the  frictional  contact  model 
[7],  is  proposed  for  the  investigation  of  mechanical  degradation 
of  CLs  under  hydrothermal  cycles,  simulating  the  start-up  and 
shutdown  of  PEM  fuel  cells. 

The  temperature  and  humidity  profiles  come  from  Ref.  [24] . 
The  driving  forces  behind  the  simulation  are  the  cycled  changes 
of  relative  humidity  (RH)  and  temperature  (7),  as  shown  in 
Fig.  12  where  the  initial  phase  difference  between  the  RH  cycle 
and  the  T  cycle  is  to  =  2  s,  and  the  period  of  humidity  and  temper¬ 
ature  are  /rh  =  tj  =  400  s.  The  maximum  and  minimum  relative 
humidity  are  RHmax  =  90%,  and  RHmin  =  30%,  respectively. 
The  maximum  and  minimum  temperature  are  Tmax  =  85  °C  and 
Tmin  =  25  °C,  respectively. 

The  2D  reconstruction  of  the  micro  structure  in  CLs  was 
shown  in  Fig.  9.  This  problem  is  also  considered  to  be  a  plane 
strain  problem  because  the  thickness  of  CLs  is  much  smaller 
than  the  planar  dimension. 

The  commercial  software  ANS  YS  is  used  to  conduct  the  sim¬ 
ulations.  The  data  used  to  simulate  the  reconstructed  CL  come 
from  the  results  of  stochastic  optimization  for  reconstruction.  At 
the  beginning  of  analysis,  the  topology  of  the  micro  structure  is 
modified  because  of  constraints  and  singularity.  Although  these 
modifications  are  very  small,  especially  in  cases  with  large  sam¬ 
ples  (i.e.,  for  large  Nx  and  Ny),  they  are  very  important  to  the 
convergence  of  computation. 

The  mesh  is  done  in  ANSYS  and  shown  in  Fig.  13.  The 
periodic  boundary  condition  is  applied  on  the  four  boundaries 
of  the  whole  computational  domain  to  avoid  the  effect  of  size. 
To  simulate  the  interface  effect  between  Nation  and  the  C/Pt 
agglomerate,  the  CZM  is  used  on  the  boundary  between  them. 
The  frictional  contact  model  is  used  to  simulate  the  contact 
phenomenon. 

Based  on  this  simulation,  Fig.  14  shows  the  evolution  of 
both  delamination  fracture  energy  and  plastic  Mises  strain. 
Fig.  14(a)  shows  two  areas  which  will  be  examined  more  closely 
in  Fig.  14(b)  and  (c)  to  investigate  the  delamination  energy  and 
plastic  strain,  respectively.  When  the  humidity  and  temperature 
begin  to  increase,  the  Nation  swells  and  delamination  fracture 
energy  accumulates.  When  Eq.  (13)  is  satisfied  as  the  result  of 
fracture  energy  accumulation  in  one  element,  this  element  is 
considered  to  be  the  failure  one.  So  as  Fig.  14(b)  shows,  the 


Fig.  12.  The  schematic  graph  of  relative  humidity  and  temperature  cycled  change. 
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Fig.  13.  Mesh  result  of  reconstructed  microstructure  with  symmetric  boundary 
condition. 

element  (indicated  with  di)  at  the  triple  boundary  of  Nation, 
C/Pt  agglomerate,  and  pore  failed  at  60  s.  Meanwhile,  the  plastic 
strain  accumulates  in  the  Nation.  At  1612  s,  there  was  crack  initi¬ 
ation  in  the  Nation  (as  indicated  with  ci  in  Fig.  14(c)).  Therefore, 
as  humidity  and  temperature  cycle,  the  plastic  strain  accumulates 
in  the  Nation,  resulting  in  fatigue  (Fig.  14(c)).  Thus  there  is 
competition  between  the  Nation’s  yield  failure  and  Nafion-C/Pt 
agglomerate  delamination  because  of  the  cyclic  environmental 
changes. 


As  we  know,  the  microstructure  based  on  stochastic  opti¬ 
mization  is  not  unique;  that  is,  a  single  microstructure  cannot  be 
determined  by  a  few  lower  order  correlation  functions,  perhaps 
because  optimization  is  not  the  true  global  minimum  and  the 
optimization  process  is  sensitive  to  the  initial  pattern  of  recon¬ 
struction.  So  different  microstructures  should  be  simulated  to 
investigate  the  common  underlying  mechanism  of  microstruc¬ 
ture  changes  in  CLs. 

However,  many  factors  in  the  micro  structure  have  complex 
influences  on  the  electrochemical  performance  of  fuel  cells. 
Phase  connectivity  is  one  of  them.  Here  phase  connectivity 
means  the  interface  between  different  solid  phases  in  CLs.  It 
is  the  key  factor  in  determining  effective  contact  area,  effective 
chemical  reaction  rate,  and  percolation  limits.  In  experiments, 
it  is  hard  to  observe  the  evolution  of  these  interfaces.  However, 
from  the  mechanical  analysis,  the  interfaces  between  Nation  and 
the  C/Pt  agglomerate  can  be  determined  by  simply  calculating 
the  contact  length  between  them. 

Five  different  microstructures  (Fig.  9)  were  simulated  under 
cyclic  changes  of  humidity  and  temperature.  The  accumulation 
of  both  delamination  energy  and  plastic  energy  were  observed 
in  these  simulations.  The  competition  between  the  accumulation 
of  delamination  energy  and  the  accumulation  of  plastic  strain 
(energy)  plays  a  key  role  in  microstructure  changes. 

As  a  result  of  the  simulation  results  for  these  five  different 
microstructures,  the  phase  connectivity  can  be  measured  after 
several  cycles,  as  shown  in  Fig.  15,  where  A L  represents  the 
change  in  the  length  of  the  interface  between  Nation  and  the 
C/Pt  agglomerate.  An  increase  in  A L  means  a  decrease  in  the 
interfaces  between  Nation  and  the  C/Pt  agglomerate.  The  con¬ 
nectivity  decreases  by  more  than  3%  of  the  whole  initial  phase 
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Fig.  14.  (a)  Elements  in  reconstructed  microstructure,  where  the  b  and  c  boxes  are  areas  which  will  be  investigated  in  subfigures  (b)  and  (c);  (b)  the  delamination 
energy  accumulation  on  the  interface  between  Nafion  and  the  C/Pt  agglomerate  at  60  s,  where  di  represents  the  initiation  position  of  the  delamination;  (c)  the  plastic 
Mises  strain  pattern  at  1612  s,  where  ci  represents  the  initiation  position  of  crack  in  the  Nafion. 
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Fig.  15.  The  evolution  of  phase  connection  length  between  Nation  and  the  C/Pt 
agglomerate. 

connectivity  after  200  cycles  for  all  five  samples.  Although  there 
are  various  sample- specific  microstructure  changes  in  the  five 
samples,  the  increase  in  connectivity  change  is  a  common  feature 
and  is  a  very  important  indication  of  the  decay  in  performance. 

These  simulations  clearly  demonstrate  that  the  mechanical 
mechanism  underlying  the  microstructure  changes  for  CLs  is 
the  competition  between  delamination  energy  accumulation  on 
the  interface  between  different  solid  phases  and  plasticity  accu¬ 
mulation  in  the  Nafion. 

5.  Conclusions 

Based  on  experimental  images  of  CLs,  statistical  features 
are  extracted  and  are  used  to  reconstruct  the  microstructure  in 
CLs.  The  reconstruction  process  is  accomplished  in  two  steps. 
One  is  the  filtering  method,  which  changes  the  Gaussian  dis¬ 
tributed  signals  to  three-phase  microstructure  reconstruction. 
The  other  is  the  stochastic  optimization  method,  which  refines 
the  microstructure  reconstruction  to  fit  the  statistical  features  of 
the  experimental  images. 

To  understand  the  underlying  mechanical  mechanism  of 
microstructure  changes,  the  reconstructed  microstructure  of  CLs 
subjected  to  hydrothermal  loading  cycles  is  investigated.  The 
mechanical  analysis  model  developed  in  our  companion  paper 
[7]  is  extended  to  calculate  the  reconstructed  microstructure 
changes.  Numerical  simulation  shows  that  delamination  and 
fracture  phenomena  do  happen  in  the  CLs  after  certain  duty 
cycles  of  humidity  and  temperature,  as  the  competition  between 
the  plasticity  energy  accumulation  in  the  Nafion  and  the  delam¬ 
ination  energy  accumulation  on  the  interface  between  Nafion 
and  the  C/Pt  agglomerate.  Although  different  reconstructed 
microstructures  show  sample  specific  results,  a  common  trend 
can  be  observed  in  the  simulation  results.  That  is,  the  cycled 


change  of  environment  causes  a  decrease  in  connection  between 
different  solid  phases,  which  may  be  an  indication  of  perfor¬ 
mance  degradation. 
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